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^ | Abstract 

The buckling of hyperelastic incompressible cylindrical shells of ar- 
\yy . bitrary length and thickness under axial load is considered within the 

\ framework of nonlinear elasticity. Analytical and numerical methods 

for bifurcation are obtained using the Stroh formalism and the exact 
solution of Wilkes [Q. J. Mech. Appl. Math. 8:88-100, 1955] for the 
linearized problem. The main focus of this paper is the range of va- 
lidity of the Euler buckling formula and its first nonlinear corrections 
qq ■ that are obtained for third-order elasticity. 

ON 
CO 

1 Introduction 

CN 

qq ' Under a large enough axial load an elastic beam will buckle. This phe- 

nomenon known as elastic buckling or Euler buckling is one of the most cel- 
ebrated instabilities of classical elasticity. The critical load for buckling was 
^ ! first derived by Euler in 1744 [HOE] an d further refined for higher modes by 

Lagrange in 1 770 jH [5] . Both authors reached their conclusion on the basis 
of simple beam equations first derived by Bernoulli [6] (See Fig. []]). Since 
then, Euler buckling has played a central role in the stability and mechanical 
properties of slender structures from nano- to macro- structures in physics, 
engineering, biochemistry, and biology [3 [8]. Explicitly, the critical compres- 
sive axial load N that will lead to a buckling instability of a hinged-hinged 
isotropic homogeneous beam of length L is 

A^Euler = ^J^, (1.1) 
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Figure 1: Euler problem: Left: illustrations from Euler manuscript [TJ. Right: 
Lagrange solutions [3], mode 1, 2, and 3. 



where "7r is the circumference of a circle whose diameter is one" [3], E is 
Young's Modulus, and I is the second moment of area which, in the case of 
a cylindrical shell of inner radius A and outer radius B, is I = 7r(5 4 — A 4 )/4. 

There are many different ways to obtain this critical value and infinite 
variations on the theme. If the beam is seen as a long slender structure, the 
one- dimensional theory of beams, elastica, or Kirchhoff rods, can be used 
successfully to capture the instability, either by bifurcation analysis, energy 
argument [7] , or directly from the exact solution, which in the case of rods can 
be written in terms of elliptic integrals [9] . The one-dimensional theory can be 
used with a variety of boundary conditions, it is particularly easy to explain 
and generalize, and it can be used for large geometric deflections of the axis 
[TO] . However, since material cross sections initially perpendicular to the axis 
remain undeformed and perpendicular to the tangent vector, no information 
on the elastic deformation around the central curve can be obtained. In 
particular, other modes of instability such as barreling cannot be obtained. 
Here, by barreling, we refer to axisymmetric deformation modes of a cylinder 
or a cylindrical shell. These modes will typically occur for sufficiently short 
structure. 

The two-dimensional theory of shells can be used when the thickness of 
the cylindrical shell is small enough. Then the stability analysis of shell 
equations such as the Donnell-von Karman equations leads to detailed infor- 
mation on symmetric instability modes, their localization and selection [11] . 
However, the theory cannot be directly applied to obtain information on the 
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buckling instability (asymmetric buckling mode). 

The three-dimensional theory of nonlinear elasticity provides, in princi- 
ple, a complete and exact description of the motion of each material point 
of a body under loads. However, due to the mathematical complexity of 
the governing equations, most problems cannot be explicitly solved. In the 
case of long slender structures under loads, the buckling instability and its 
asymptotic limit to Euler criterion can be captured by assuming that the 
object is either a rectangular beam [121 H21 H3j or a cylindrical shell under 
axial load. Then, using the theory of incremental deformations around a 
large deformation stressed state, the buckling instability can be recovered by 
a bifurcation argument, usually referred to, in the nonlinear elasticity theory, 
as small- on-large, or incremental, deformations. In comparison to the one- 
and two-dimensional theories, this computation is rather cumbersome as it 
is based on non-trivial tensorial calculations, but it contains much informa- 
tion about the instability and the unstable modes selected in the bifurcation 
process. 

Here, we are concerned with the case of a cylindrical shell under axial load. 
This problem was first addressed in the framework of nonlinear elasticity in a 
remarkable 1955 paper by Wilkes [15] who showed that the linearized system 
around a finite axial strain can be solved exactly in terms of Bessel functions. 
While Wilkes only analyzed the first axisymmetric mode (n = 0, see below), 
he noted in his conclusion that the asymmetric mode (n = 1) corresponds 
to the Euler strut and doing so, opened the door to further investigation 
by Fosdick and Shield [16] who recovered Euler's criterion asymptotically 
from Wilkes's solution. These initial results constitute the basis for much 
of the modern theory of elastic stability of cylinders within the framework 
of three-dimensional nonlinear elasticity [171 EIH1 CHI EDI EI]- The experimen- 
tal verification of Euler's criterion was considered by Southwell [22J and by 
Beatty and Hook |23j . 

The purpose of this paper is threefold. First, we revisit the problem of 
the stability of an incompressible cylindrical shell under axial load using the 
Stroh formalism [21] and, based on Wilkes solution, we derive a new and 
compact formulation of the bifurcation criterion that can be used efficiently 
for numerical approximation of the bifurcation curves for all modes. Second, 
we use this formulation to obtain nonlinear corrections of Euler's criterion 
for arbitrary shell thickness and third-order elasticity. Third, we consider the 
problem of determining the critical aspect ratio where there is a transition 
between buckling and barreling. 
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2 Large deformation 



We consider a hyperelastic homogeneous incompressible cylindrical tube of 
initial inner radius A, outer radius B, and length L, subjected to a uniaxial 
constant strain A3, and deformed into a shorter tube with current dimensions 
a, b, and I. The deformation bringing a point at (R, B, Z), in cylindrical co- 
ordinates in the initial configuration, to (r, 9, z) in the current configuration 
is 

r = A1.R, 9 = 9, 2 = X 3 Z, (2.1) 

where Ai = a/ A = b/B and A 3 = l/L. The physical components of the 
corresponding deformation gradient F are 

[F] = diag(A 1 ,A 1 ,A 3 ), (2.2) 

showing that the principal stretches are the constants Ai, A 2 = Ai, A 3 ; and 
that the pre-strain is homogeneous. Because of incompressibility, det F = 1, 
so that 

Ai = A 3 1/2 . (2.3) 
The principal Cauchy stresses required to maintain the pre-strain are [25] 

dW 

ct^-p + A;— , (i = 1,2,3), (2.4) 

(no sum) where p is a Lagrange multiplier introduced by the internal con- 
straint of incompressibility and W is the strain energy density (a symmetric 
function of the principal stretches). In our case, a 2 = <J\ because A2 = Ai. 
Also, o\ = because the inner and outer faces of the tube are free of traction. 
It follows that 

p = XtWx, a 3 = X 3 W 3 - A1W1, (2.5) 

where Wi = dW/dXi, and we conclude that the principal Cauchy stresses are 
constant. 



3 Instability 

To perform a bifurcation analysis, we take the view that the existence of small 
deformation solutions in the neighborhood of the large pre-strain signals the 
onset of instability [26J. 
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3.1 Governing equations 

The incremental equations of equilibrium and incompressibility read [25] 

div s = 0, div u = 0, (3.1) 

where s is the incremental nominal stress tensor and u is the infinitesimal 
mechanical displacement. They are linked by 

s = A® grad(w) + p grad(w) — pi, (3.2) 

where p is the increment in the Lagrange multiplier p and ,4o is the fourth- 
order tensor of instantaneous elastic moduli. This tensor is similar to the 
stiffness tensor of linear anisotropic elasticity, with the differences that it 
possesses only the minor symmetries, not the major ones, and that it reflects 
strain-induced anisotropy instead of intrinsic anisotropy. Its explicit non-zero 
components in a coordinate system aligned with the principal axes are [25J: 

•Aoiijj Aj Xj ^Vij , 

A^(W i -A J ^)A i 2 /(A 4 2 -A J 2 ) 

<Aoijij\>Aoiiii ^Oiijj A«Wj)/2 
•Aoijji >Aojiij>Aoijij AjPVj, 

(no sums) where Wij = d 2 W / {dXidXj) . Note that some of these components 
are not independent because here Ai = A2. In particular we have 

.A02121 = .4.01212, .4.02323 = .4.01313, *4o2222 = -4omi, (3-4) 

>4o2233 = >4oil335 -4o2332 = -4oi331; *4o3232 = -4o3131) 

*4oi221 + -4oi212-4oilll — .4.01122 = 2^4oi212 + «4oi331 — *4oi313- 

3.2 Solutions 

We look for solutions which are periodic along the circumferential and axial 
directions, and have yet unknown variations through the thickness of the 
tube, so that our ansatz is 

{u,v,w,p, s rr ,s r9 ,s rz } = 

{U(r), V(r), W(r),P(r),S rr (r), S rd (r), S rz (r)}e^ ne+kz \ (3.5) 



if i ^ j, Xi ^ Xj, 
if i ^ j ; Ai = Xj, 



(3.3) 
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where n = 0, 1, 2, . . . is the circumferential number, k is the axial "wave"- 
number, and all upper-case functions are functions of r alone. 

The specialization of the governing equations (13.11) to this type of so- 
lutions has already been conducted in several articles (see for instance [151 
EES1 EES [23 [2H] , and Dorfmann and Haugton [21] for the compressible coun- 
terpart.) Here we adapt the work of Shuvalov [2H] on waves in anisotropic 
cylinders to develop a Stroh-like formulation of the problem [21]. The central 
idea is to introduce a displacement-traction vector, 



>rz\ i 



rj = [U, V, W, irS rr , irS r g, irS r 
so that the incremental equations can be written in the form 

-^-r/(r) = -G(r)ri(r), 
dr r 

where G is a 6 x 6 matrix, with the block structure 



G 



G\ (?2 



G 2 — G 2 , 



Gz 



(3.6) 



(3.7) 



(3.8) 



Here the superscript '+' denotes the Hermitian adjoint (transpose of the 
complex conjugate) and G\, G 2 , G 3 are the 3x3 matrices 



i — n 
—n — i 

— kr 



-kr 










AH212 










A 



respectively, with 

Kll = 4.4.01212 + (-4.03131 
K 12 = 4^^.01212, ^13 
«33 



- A 01313 )k 2 r 2 , 



2A 



01212 



kr. 



01313- 

«22 = 



Kii IK12 Kl3 
-IK.12 K,22 ^23 
-iKl3 K23 K,33_ 

(3.9) 

= n(2.4 i2i2 + A 01 3 13 )kr, 

4n 2 Al212 + .4.03131 



n2 ^01313 + (4^4.01212 + 2^4.01122 — -4oilll + -4-03333 — 2^4oil33)^ 2r2 - 



(3.10) 



As it happens, there exists a set of 6 explicit Bessel-type solutions to 
these equations when n ^ 0. This situation is in marked contrast with 
the corresponding set-up in linear anisotropic elastodynamics, where explicit 
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Bessel-type solutions exist only for transversely isotropic cylinders with a 
set of 4 linearly independent modes, and do not exist for cylinders of lesser 
symmetry (301 EH]- As mentioned in the Introduction, the 6 Bessel solutions 
are presented in the article by Wilkes [15] (for a derivation see Bigoni and 
Gei [20]). 

First, denote by qf, q 2 the roots of the following quadratic in q 2 , 

-4.01313? 4 — (-4-01111 + -4-03333 — 2.4.01331 — 2*4.oii 33 )g 2 + -4.03131 = 0. (3.11) 

Then the roots of this quadratic are ±gi and ±52, and it can be checked that 
the following two vectors are solutions to (13. 7p . 



^(1)^(2) 



Tl 

il' n {qkr), —l n (qkr), -ql n (qkr) 

qkr 



(Ai3i3<? 2 + A)^i)l n {qkr) + 2A1212 In(# r ) ^-Iniqkr) , 

q \ qkr J 

-2mAi2i2 (l' n (qkr) - -^l n (qkr)^j , -i(l + q^krA^zMqkr) , (3.12) 

where q = qi,q 2 m turn and I n is the modified Bessel function of order n. 
Similarly, we checked that the following vector tj^ 3 \ 



71 

ii-in(g3^), -qM^hr), 0, 

L kr 



2nq 3 A 012 i2 (l' n {qzhr) l —l n {q^kr) ] 

V q3kr J 



ig 3 Ai2i2 (2I' n (q 3 kr) - (q 3 kr + 2 



n 2 



l n (q3kr) 



q 3 kr y 

-in^4oi3i3l n (g 3 A;r)]* , (3.13) 

is also a solution, where q 3 is root to the quadratic 

Ai2i2<? 2 - Asm = 0. (3.14) 

Finally we also checked that the vectors r)^ A \ r}^ 5 \ and 77^, obtained by 
replacing I n with the modified Bessel function K n in the expressions above, 
are solutions too. 

Next, we follow Shuvalov [29] and introduce Af(r) as a fundamental ma- 
trix solution to (13.71) . 

M{r) = [t7 (1) |77 (2) |...|t7 (6) ] . (3.15) 



It clearly satisfies 

T-A/*(r) = -G(r)Af(r). (3.16) 
dr r 

Let M(r, a) be the matricant solution to (13. 7ft . that is the matrix such that 

r/(r) = M(r, 0)17(0), M(a, a) = I( 6 )- (3.17) 

It is obtained from A/"(r) (or from any other fundamental matrix made of 
linearly independent combinations of the r]^) by 

M(r, a) = Af^Af-^a), (3.18) 

and it has the following block structure 

Mx(r,a) M 2 {r,a)~ 
M 3 (r,a) M 4 (r,a)_ 



M(r, a 
say. 



(3.19) 



3.3 Boundary conditions 

Some boundary conditions must be enforced on the top and bottom faces of 
the tubes. Considering that they remain plane (W = on z = 0, /) and free 
of incremental tractions (S rz = S r g = on z = 0, 1) leads to 

mTi mix 

where m = 1, 2, 3, . . . but, since the equations depend only on k, we can take 
m = 1 without loss of generality. 

The other boundary conditions are that the inner and outer faces of the 
tube remain free of incremental tractions. We call S = [S rr , S r g, SVz]* the 
traction vector, and U = [U, V, W]* the displacement vector. We substitute 
the condition S(a) = into (13. IT)) and (13. 19)) to find the following connection, 

rS(r) = z(r,a)U(r), where z = iM 3 Mf 1 (3.21) 

is the (Hermitian) 3x3 impedance [29J. Since S(b) = 0, a non-trivial solution 
only exists if the matrix z(b,a) is singular, which implies the bifurcation 
condition 

d6t ZM - 1 detM 1 (b,a) =°- (3 ' 22) 

This is a real equation because z = z + [29]; note that the nature (i.e. real 
or complex [TjJ], simple or double [2T]) of the roots gi, qz is irrelevant. 
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4 The adjugate method 



We are now in a position to use the bifurcation condition (13.221) to compute 
explicitly bifurcation curves for each mode n. We note that the components 
of ^4o depend on the strain energy density W and on the pre-strain, which by 
(12.31) depends only on A3; so do <Zi, ?2> (73 by (13.111) and (13.141) . According to 
(I3.12p and (13.131) . the entries of M(b, a) thus depend (for a given W) on A3, 
n, ka, and kb only. For a given material (W specified) with a given thickness 
(b/a = B/A specified), the bifurcation equation (13.221) gives a relationship 
between a measure of the critical pre-stretch: X^X^ 2 , and a measure of the 
tube slenderness: kb = 2-RmibjV) = 27cm\ 3 3 ^ 2 (B/L), for a given bifurcation 
mode (n specified). That is, for a given tube slenderness, what is the axial 
strain necessary to excite a given mode. 

While this bifurcation condition is formally clear, it has not been success- 
fully implemented to compute all bifurcation curves. Indeed, for mode n > 1, 
the numerical root finding of det(z) becomes unstable (as observed in [21] for 
a similar problem) and, in explicit computations, most authors do not use the 
exact solution by Wilkes but use a variety of numerical techniques to solve 
the linear boundary value problems directly (such as the compound matrix 
method [32J, the determinantal method [33], or the Adams- Moulton method 
|34|). Note that from a computational perspective, the Stroh formalism is 
particularly well-suited and well-behaved [35j [36] and if numerical integra- 
tion was required it would provide an ideal representation of the governing 
equation. 

Rather than integrating the original linear problem numerically, we now 
show how to use an alternative form of (13.221) to compute all possible bi- 
furcation curves. This method bypasses the need for numerical integration 
and reduces the problem to a form that is manageable both numerically and 
symbolically, to study analytically particular asymptotic limits. The main 
idea is to transform condition (13.221) by factoring non-vanishing factors. We 
start by realizing that since the fundamental solutions {rj^^i = 1 . . . , 6} are 
linearly independent, the matrix Af(r) is invertible for all r e [a, b], which 
implies that the elements of M(r,a) are bounded for r e [a, b}. Therefore 
det(Mi(r,a)) is bounded and detz = implies det(M" 3 (6, a)) = 0. Instead 
of expressing det(M 3 (b, a)) as the determinant of a 3 x 3 submatrix of a 
matrix obtained as the product of two 6x6 matrices, we first decompose 
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A/"(r) as 

^ |M(r) AA 2 (r) 



(4.1) 



say, where each block is a 3 x 3 matrix. We also rewrite Eq. (13.181) as 

M(r, a)Af(a) = Af(r). (4.2) 

and write explicitly the two entries A/" 3 (r) and A/" 4 (r), which are 

M 3 (r, a)AA(a) + M 4 (r, a)Af 3 (a) = A/" 3 (r), (4.3) 

M 3 (r, a)Af 2 (a) + M 4 (r, a)Af 4 (a) = A/" 4 (r), (4.4) 

which implies 



M 3 (r,a) = 

[A^A^a) " A^A^a)] [JSfiiaW^W-rffaWl 1 * 



a 



-1 

(4.5) 



Using again the fact that the entries of J\f are bounded, we have that the 
bifurcation condition det(A / t 3 (6, a)) = implies that 

detQ(6,a) = (4.6) 

where 

Q(b,a) = det(A/>))A^)adj(Af 3 (a)) - det(Af 3 (a))Af 4 (&)adj(A^ 4 (a)), 

(4.7) 

and adj(A) is the adjugate matrix of A, that is the transpose of the cofactor 
matrix (which in the case of an invertible matrix is simply adj(A)=det(A) A -1 ). 
This new bifurcation condition is equivalent to the previous one but has many 
advantages. The matrix Q only involves products of 3 x 3 matrices and is 
polynomial in the entries of Af, that is, det Q(b, a) is a polynomial of degree 
18 in Bessel functions and has no denominator (hence no small denomina- 
tor). Both numerically and symbolically, this determinant is well-behaved, 
even in the limits a — > 0, which corresponds to a solid cylinder, and n = 0, 
which corresponds to the first barreling mode (and usually require a special 
treatment). We will refer to the use of this form of the bifurcation condition 
as the adjugate method. 
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4.1 Numerical results 



As a first test of the stability of the numerical procedure, we consider a neo- 
Hookean potential W = C\(Ii — 3)/2, for the typical values C\ = 1, B/A = 2 
and compute (and plot in Fig. [2]) for the ten first modes (n = to n = 9), 
the critical value of A = A3 as a function of the current stubbiness kb = irb/l 
(the initial stubbiness is v — B/L = kb\ 3 ^ 2 /ir). The known classical fea- 
tures of the stability problem for the cylindrical shell are recovered, namely: 
for slender tubes, the Euler buckling (n = 1) is dominant and becomes un- 
avoidable as the slenderness decreases; there is a critical slenderness value at 
which the first barreling mode n = is the first unstable mode (in a thought 
experiment where the axial strain would be incrementally increased until the 
tube becomes unstable); and for very large kb, the critical compression ratio 
tends asymptotically to the value A = 0.444, which corresponds to surface 
instability of a compressed half-space [12] . 




Figure 2: Bifurcation curves (stretch as a function of stubbiness) of an ho- 
mogeneous neo-Hookean cylindrical tube for modes n = to n = 9 with 
b/a = B/A = 2 and C x = I. 

For a second test, we consider very thin neo-Hookean tubes with B/A = 
1.01. Here we are interested in the mode selection process. As the stubbiness 
increases, the buckling mode rapidly ceases to be the first excited mode and 
is replaced by different barreling modes. From Fig. [3j it appears clearly 
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that as kb increases, modes n = 1 to 9 are selected (modes n = and 
n = 10 remain unobservable) . There is one particularly interesting feature in 
these two sets of bifurcation curves. Depending on both the tube thickness 
and stubbiness, the instability mode of a tube transition from buckling to 
barreling, the material transition from either the one-dimensional behavior of 
slender column to the two-dimensional behavior of a thin short tube, or the 
three-dimensional behavior of a thick short tube. Accordingly we will refer 
to these particular geometric values where transition occurs as dimensional 
transitions and obtain analytical estimates for them in the next section. 
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Figure 3: Bifurcation curves (stretch as a function of stubbiness) of an ho- 
mogeneous neo-Hookean cylindrical tube for modes n = to n = 10 with 
b/a = B/A = 1.01 and C x = 1. 



5 Asymptotic Euler buckling 

We are now in a position to look at the asymmetric buckling mode (n = 1) 
corresponding to Euler buckling in the limit A — > 1. The asymptotic form 
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of Euler criterion cannot be obtained for a general strain-energy density. 
This is why we choose the Mooney-Rivlin potential, which for A close to 1, 
corresponds to the most general form of third-order incompressible elasticity 
(see next Section), 

W = Ci(Ji-3)/2 + C 2 (/ 2 -3)/2, (5.1) 

where C\ > and C 2 > are material constants, I\ = Af + A| + A§, and 
I 2 = A^A 2 , + AgA| + A|Af. Close to A = 1, we introduce a small parameter 
related to the stubbiness ratio 

e = kb = Tib/l, (5.2) 
and look for the critical buckling stretch A as a function of e of order 2M: 

M 

A = A(e) = 1 + J2 A ^ e2m + 0(^ 2M+2 )- (5-3) 

m=l 

Similarly, we expand d(X) = det Q(b, a) in powers of e, 

M d 

d(\) = J^d 2m e 2m + 0(e 2M ^ 2 ), (5.4) 

m=l 

and solve each order d 2m = for the coefficients A 2m . This is a rather cum- 
bersome computation. The first non-identically vanishing coefficient appears 
at order 24 and a computation to order 28 is necessary to compute the correct 
expression for A, which is found to be to order 6 in e 



A = 1 + A (2) e 2 + A (4) e 4 + A (6) e 6 + 0(e 8 ), (5.5) 



with 



a (2 ) = -4^> ( 5 - e ) 



(19C 2 + 28C 1 )p 4 + 2 (53C 2 + 62C 1 )p 2 + 19C 2 + 28C 1 
A(4) " 144(d + C 2 )p 4 ' (5 ' 7) 

A(6) = ~A608p^-l)(C 1 + C 2 ) X (5 ' 8) 
(973 d + 341 C 2 ) (p 10 - 1) + (7073 d + 3385 C 2 ) p 2 (p 6 - 1) + 

43921n(p)(C 1 + C 2 )(p 6 + p 4 ) + 4(377C 2 + 1141C 1 )(p 6 - p 4 
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where p = B j A = b/a. It is of interest to compare the different approxima- 
tions. We recover Euler formula by keeping only the term up to e 2 which we 
denote by Euler2. We define similarly Euler4 and Euler 6 by keeping terms 
up to order 4 and 6 in e. In FigJU we show the different approximations as 
a function of e 2 for p = 1.01 (on the left) and for p = 10 (on the right). The 
classical Euler formula is well-recovered in the limit e — > but the Euler 4 and 
Euler 6 approximations clearly improve the classical formula for larger values 
of e. It also appears from the analysis of Euler4 that for Ci ^ the classical 
Euler formula always underestimates the critical stretch for instability. 




Figure 4: Comparison of the different Euler formulae obtained by expanding 
the exact solution to order 2 (classical Euler buckling formula), Euler 4 and 
Euler 6 for a neo-Hookean potential C\ — 1, C% — 0. For comparison purpose, 
we show the critical stretch for mode n = 1 versus e 2 in which case the graph 
becomes linear in the limit e — > 0. Left: p = b/a = B/A = 1.01, Right: 
p = b/a = B/A = 10. 



6 Nonlinear Euler buckling for third-order elas- 
ticity 

The analytical result presented in the previous section was formulated in 
terms of parameters and quantities natural for the computation and the 
theory of nonlinear elasticity. In order to relate this result to the classical 
form of Euler buckling, we need to express Eq. (I5.5P in terms of the initial 
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geometric values A, B, L, the axial load acting on the cylinder, and the elastic 
parameters entering in the theory of linear elasticity. 

We first consider the geometric parameters. We wish to express the crit- 
ical load as a function of the initial stubbiness v = B/L and tube relative 
thickness p = B/A. Recalling that e = nb/l and that A = l/L,b = \~ l l 2 B, 
we have 

e 2 A 3 = vrV. (6.1) 

To express e as a function of z/, we expand e in powers of v to order 6, and 
solve (16. ip to obtain 

e 2 = vrV - 37T 4 A (2) z/ 4 - (3vr 6 A (4 ) - 15tt 6 A 2 2) )z/ 6 , (6.2) 

where A( 2 ) and A( 4 ) are defined in (I5.6H5.7I) and come from the expansion of 
A in powers of e. 

Second, we want to relate the axial compression to the actual axial load 
N. To do so, we integrate the axial stress over the faces of the tubes that is 

N = -2tt f ra 3 dr. (6.3) 

J a 

Since a 3 is constant and given by (12.51) we have 

N = -n(b 2 - a 2 )a 3 = -^(B 2 - A 2 )a 3 

A 

= - J(5 2 - A 2 ) [(A 4 - A)Cx + (A 3 - 1)C 2 ] . (6.4) 

Third, we relate the elastic Mooney parameters C\ and C 2 to the classical 
elastic parameters. Here, we follow Hamilton et al. [371 [38] and write the 
strain-energy density to third-order for an incompressible elastic material as 

W = -2{ii 2 + n 3 i 3} (6.5) 

where /i is the usual shear modulus, or second Lame parameter, and n 3 is a 
third-order elasticity constant; /i is related to Young's modulus by E = 3/i; 
also, in Murnaghan's notation, n 3 = n and in Landau's notation, n 3 = 
A, see Norris [39] for other notations. In (16.51) . ii,i 2 ,i3 are the first three 
principal invariants of the Green-Lagrange strain tensor, related to the first 
three principal invariants Ii, J 2 , /3 of the Cauchy-Green strain tensor by 

h = 2ii + 3, I 2 = 4ix + Ai 2 + 3, I 3 = 2ii + 4i 2 + 8i 3 + 1, (6.6) 
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Since 7 3 = 1, we can solve this linear system for i 2 and i 3 and write the 
strain-energy density (16. 5p as a function of l\ and J 2 , that is 



which by comparison with ( 15. ip leads to 

C 1 = 2 f i + 7i 3 /4, C 2 = [i- n 3 /4. (6.8) 

To write the nonlinear buckling formula, we consider (16 .4p and first ex- 
pand A in e using (I5.5p . then expand e in z/ using (16.21) . and, finally, substitute 
the values of the moduli in terms of the elastic parameters, which yields 

3^V(/-l)^ 

4 p 4 v ; 

1 n 5 B 2 (p 2 - 1) (20 p 4 p + 9 p 4 n 3 + 176 p 2 p + 18 p 2 n 3 + 20 p + 9 n 3 ) u A 
96 p 1 + 

5W 



512p 8 p(p 2 + 1) 

323 p 2 p 8 - 3 n 3 2 - 240 p 2 p n 3 - 9 p 2 n 3 2 - 9 p 10 p 2 + 9 p 2 + 3 p 10 n 3 2 + 
1464 In (p) p 6 p 2 + 1464 In (p) p 4 p 2 + 240 p 8 p n 3 + 240 p 6 p n 3 - 
180 p 6 p 2 + 180 p 4 p 2 + 9 p 8 n 3 2 + 6 p 6 n 3 2 - 6 p 4 n 3 2 - 323 p 2 p 2 - 240 p 4 p n 3 

While it is not surprising, it is comforting to recover to order v 2 the classical 
Euler buckling formula (11.11) (using p = B/A, v = B/L, and p = E/3). 



7 Dimensional transition 

Finally, we use the buckling formula to compute the transition between modes 
as parameters are varied. That is, to identify both the geometric values and 
the axial strain for which there is a transition between buckling and barreling 
modes. Here we restrict again our attention to the neo-Hookean case (with 
C\ = 1). From Fig. [2] and Fig. [31 it appears clearly that for e small enough 
there is a transition (depending on the value of p) from either mode n = 1 
to mode n = (large p), or from mode n = 1 to mode n = 2 (p close to 
1) as e increases. We refer to this transition as a dimensional transition, 
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in the sense that the material mostly behaves as a slender one-dimensional 
structure when it buckles according to mode n = and mostly as a two- 
dimensional structure when it barrels with mode n = 2. Indeed both modes 
of instability can be captured by, respectively, a one- or a two-dimensional 
theory. For p close to unity, the transition n = — > n = 1 occurs for small 
values of e. Therefore in this regime, we can use the approximation (15.51) 
for the barreling curve and substitute it in the bifurcation condition of mode 
n = 2. Expanding again this bifurcation condition in e as well as p, one 
identifies the values pt of p and At of A at which the transition occurs, as 

pt = i + 3 e 2_ 53 £ 4 2393 = x _l 2 13 4 _665 6 + 0(e8) 

H 4 32 384 y h 4 8 96 V ; 

(7.1) 

In terms of the initial stubbiness v = B/L, we have 



p t = l + - 7T V - - vrV + — ttV 6 + O (is 8 ) (7.2) 
Ht 4 32 384 v > K ' 

This relationship also provides a domain of validity for the Euler buckling 
formula. For sufficiently slender tube (y small), the buckling mode disappears 
when p > p t at the expense of the n = 2 barreling mode. For stubbier 
and fuller tubes, this approximation cannot be used. To understand the 
dimensional transition, we solve numerically the bifurcation condition, using 
the adjugate method, for the intersection of two different modes. That is, 
for a given value of p*, we find the value of e* such that both the bifurcation 
for either modes n = 1 and n = 2, or modes n = 1 and n = are satisfied. If 
the corresponding value A* is the largest value for which a bifurcation takes 
place, the pair (e^p*) is a transition point. The corresponding transition 
point in terms of the initial parameters is (u* = — (A*) 3 / 2 , p*). In Fig. [5] we 
show a diagram of all such pairs for both transitions. 



8 Conclusion 

This paper establishes a reliable and effective method to study the stability 
of tubes based on the exact solution of the incremental equations proposed 
by Wilkes [15] within the Stroh formalism. It then puts the method to 
use, to obtain the first geometric and material corrections to Euler buckling. 
The method can be also used to obtain the transition between buckling and 
barreling modes when a tube becomes unstable. 
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Figure 5: Dimensional transition for a neo-Hookean cylindrical tube of initial 
length L and initial radii A and B. All tubes in the n — 1 regions will become 
unstable by buckling. As the tubes get stubbier or thinner (arrows), it will 
not buckle but instead will be subject to a barreling instability Note that 
only the transition curves from mode n = are shown. Tubes in the barreling 
regions may be subject to other unstable modes. 

The method presented here can be easily generalized to different materials 
and different boundary conditions. For instance, using the exact solution of 
the incremental equations proposed in [21] for compressible materials and 
the adjugate method, an explicit form of the bifurcation condition in terms 
of Bessel functions can be obtained by following the steps presented here 
and various asymptotic behaviors can be obtained. Similarly, a variety of 
boundary value problems can be analyzed by the adjugate method, such 
as the stability problem of a tube under pressure and tension (3H HD] , the 
problem of a tube embedded in an infinite domain [2U], and the problem of 
a tube with coating [H]. In all these cases, useful asymptotic formulae for 
the buckling behavior could be obtained by perturbation expansions. 

It is also enticing to consider the possibility of performing an analytical 
post-buckling analysis of the solutions. Since the solutions of the linearized 
problem can be solved exactly, a weakly non-linear analysis of the solution 
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should be possible to third-order. This would yield, in principle, an equation 
for the amplitude of the unstable modes containing much information about 
the actual amplitude of the unstable modes but also on the localization of 
unstable modes after bifurcation. We leave this daunting task for another 
day. 
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